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A new meron cluster algorithm is constructed to study the finite temperature critical behavior of the chiral condensate in a (3 + 1) 
dimensional model of interacting staggered fermions. Using finite size scaling analysis the infinite volume condensate is shown to 
be consistent with the behavior of the form (T c — y) 314 ( 7 ) f or temperatures less than the critical temperature and m 1 ^ 4 ' 87 ' 10 ^ at 
the critical temperature confirming that the critical behavior belongs to the 3-d Ising universality class within one to two sigma 
deviation. The new method, along with improvements in the implementation of the algorithm, allows the determination of the 
critical temperature T c more accurately than was possible in a previous study. 



1. Motivation 

The construction of Monte Carlo algorithms to solve 
problems in many body quantum mechanics involving 
fermions is notoriously difficult. This difficulty is reflected 
in our inability to perform precise quantitative calcula- 
tions in strongly interacting fermionic models which are 
necessary to understand a variety of phenomena includ- 
ing high temperature superconductivity and the physics of 
strongly interacting matter. The essential problem arises 
due to the Pauli principle which can produce negative 
Boltzmann weights when the quantum partition function 
is rewritten as a path integral in a convenient basis. As 
a result, the probability distribution that should be used 
for importance sampling is unclear. 

The conventional approach to problems involving 
fermions is to integrate them out in favor of a fermion 
determinant. In cases where this determinant is positive 
it is often possible to use known sampling methods for a 
bosonic problem to devise an algorithm |lj-f|]. Some of 
these methods are inexact since they involve discretiza- 
tion of a differential equation and require some care and 
study before they can be applied to a new problem. Oth- 
ers can suffer from truncation errors that approximate the 
original partition function. Worst of all, these algorithms 
suffer from the usual problems of critical slowing which 
makes it difficult to study phase transitions using them. 

The study of phase transitions, especially in the context 
of fermionic models, is of interest in a variety of fields. 
In condensed matter physics strong correlations between 
electrons can lead to many interesting critical effects like 
high T c superconductivity @] and quantum phase transi- 
tions . In high energy physics, existence of new phases 
and fixed points have been predicted which may lead 
to novel formulations of quantum field theories beyond 



perturbation theory. Additionally, exotic phases arise in 
dense nuclear matter due to strong interactions among 
quarks |ic|| . 

Since fermions acquire a screening mass on the order 
of the temperature T, one expects that the finite tem- 
perature critical behavior close to a second order phase 
transition in a (d + 1) dimensional theory is governed by 
a d dimensional low energy effective theory that is purely 
bosonic [p"T|| . A few years ago, this conventional wisdom 
was questioned based on a large N calculation jl^] in a 
Gross-Neveu model. It was shown that the finite temper- 
ature phase transition in the (2 + 1) dimensional model re- 
produced mean field exponents instead of the expected 2-d 
Ising exponents. This claim was later backed by numeri- 
cal evidence |l3| using the hybrid Monte Carlo algorithm. 
Although, the reason for the unexpected critical behavior 
was later attributed to the narrowness of the Ginsburg 
region in the large N limit the numerical evidence 
provided to substantiate the earlier claims remains dis- 
turbing and perhaps shows the inadequacy of the numer- 
ical methods used. Recently, the chiral transition in two 
flavor QCD with an additional four-fermion interaction 
was studied using the hybrid molecular dynamics algo- 
rithm which showed evidence for non-mean field critical 
exponents | fl5"|| . However, again the expectations based on 
simple dimensional reduction were not observed, instead 
the data appeared to be consistent with tricritical behav- 
ior. 

The inability to provide conclusive answers to ques- 
tions related to the critical behavior in fermionic theo- 
ries is closely related to the lack of efficient fermion algo- 
rithms. The fermion cluster algorithms, recently proposed 
in |l6|,[l^] , provide a novel approach to the problem. The 
new method, referred to as the meron cluster algorithm, 
uses well known quantum cluster algorithm techniques Jig] 
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to solve the fermion sign problem completely. The first ap- 
plications of the meron algorithm emerged last year when 
it was used to study the critical behavior in a relativistic 
system of interacting staggered fermions with a discrete 
chiral symmetry [^9|2C| l. The results indicated that the 
symmetry is broken by the ground state but is restored 
by thermal fluctuations at high temperatures. However, 
in order to avoid the complications that arise in the al- 
gorithm due to the addition of a mass term, the previ- 
ous study focused on massless fermions. Since the chiral 
condensate vanishes in this case the scaling of the chiral 
susceptibility with the volume was used to find the critical 
temperature and the critical exponents v and 7. 

In this article a new meron algorithm is proposed and 
applied to study the staggered fermion model in the pres- 
ence of the mass term[]. This makes it possible to study the 
critical behavior of the chiral condensate. The main result 
of the article is that in (3+1) dimensions the chiral con- 
densate behaves like (ijrtp) = A(T C — T)^ just below T c with 
(3 — 0.314(7) and vanishes for higher temperatures. The 
corresponding exponent in the 3-d Ising model is known 
to be 0.32648(18). As a bonus the critical temperature 
is also determined more accurately than previous stud- 
ies. Furthermore at T c the mass dependence of the con- 
densate also behaves as expected, (pjnp) = Bm 1 / 6 , where 
S = 4.87(10) as compared to 4.7893(22) of the Ising model. 
Conventional fermionic algorithms have never been able to 
confirm the predictions of universality in a strongly inter- 
acting fermionic model with such precision. 

2. The Model 

The Hamilton operator for staggered fermions hopping 
on a 3-d cubic spatial lattice with V = L 3 sites (L even) 
and anti-periodic spatial boundary conditions, considered 
here, is given by 

H = J~]h Xt i +my^s I . (1) 

x : i x 

The term 

h x ,i = ^(4c x+i + c\ + f x ) + {n x - i)(n x+ j - ~), (2) 

couples the fermion operators at the nearest neighbor sites 
x and x + i where i is a unit-vector in the positive i- 
direction and the mass term 

s x = (-l)* 1+ * a+ * a (n t -l/2), (3) 

is a single site operator. The fermion creation and an- 
nihilation operators c x and c x satisfy the canonical anti- 
commutation relations and n x = c x c x is the number op- 
erator. The phase factors T) Xt i = 1, rj X; 2 = (— l) Xl and 

1 Preliminary results of this work was presented in 1221 



Vx,3 = (— l) Xl+X2 are well known in the staggered fermion 
formulation f2lf . 

This model was originally studied in the chiral limit 
(rn = 0) in jis|] . In this limit the Hamilton operator 
is invariant under shifts in the X3 direction. The mass 
term breaks this symmetry up to shifts by an even num- 
ber of lattice units, thus breaking a Z2 symmetry which 
can be related to a subgroup of the well known chiral 
symmetry of relativistic massless fermions. This symme- 
try is broken spontaneously at zero temperatures, while 
thermal fluctuations restore it at some high temperatures 
|ljj. In this article the critical behavior near the sec- 
ond order transition is studied using the chiral condensate 
(ipip) = — J2x( s x)/V usm g the algorithm presented in p2( | 
and briefly sketched here. 

The construction of the path integral for the parti- 
tion function is standard. First, the Hamilton operator 
is rewritten as H = Hi + H2 + ■ ■ ■ + Hq + with 

Hi = ^ h x ,i H i+3 = ^ h x ,i (4) 

Xi even Xi odd 

for i = 1,2,3 and H-j — rn^2 x s x . The partition function 
is then approximated by 

Tr[e~ eHl e~ eHr ^ 6 e~ eH2 e~ eH7 ^ 6 .. _ e - eH e e -^H 7 /6-iM ^ 

where the inverse temperature has been divided into M 
slices such that Me = 1/T. At a fixed temperature, the 
above approximation becomes exact in the limit M — > 00 
and e — > 0. On the other hand for any fixed M, the 
approximation defines a new theory with a phase structure 
and critical behavior that can be identical to the M = 00 
theory. For simplicity this article focuses on the theory 
with M — 4. 

3. Meron Cluster Algorithm 

In order to solve the model discussed in the previous sec- 
tion using a meron cluster algorithm, the partition func- 
tion should first be written in terms of fermion occupation 
numbers n and bond variables b so that it can be repre- 
sented by 

Z = Siga[n,b]W[n,b]. (6) 

[n,5] 

Configurations [n, b] live on an L d x 4dM lattice. A typ- 
ical configuration in one space and one time dimension is 
shown in figure [l| The shaded regions represent either 
the nearest neighbor interactions that arise due to terms 
of the form exp(— th x ^) or the single site interactions due 
to mass terms exp(— ems x /6). The various steps that are 
used to determine W[n, b] and Sign[n, b] starting from the 
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Figure 1 . A configuration of fermion occupation numbers 
and bonds on a one dimensional periodic spatial lattice of 
size L = 4. Each shaded region corresponds to a transfer 
matrix clement and has weights and signs that are shown 
in figure ||. 



Figure 2. The magnitude and sign of the non-zero trans- 
fer matrix elements. The product of such weights over 
the entire lattice determine W[n, b] and Sign[n, b}. The 
signs £ get contributions from local staggered fermion 
phase factors and non-local factors that arise due to anti- 
commutation relations. 



partition function (||) have been discussed in ]2l^ , |l9| , p4j . 
The final results can be represented as a set of simple 
rules. W[n, b] turns out to be a product of magnitudes 
of the transfer matrix elements and Sign[n, b] represents 
the product of their signs. The non-zero matrix elements 
are shown in figure ^ along with their weights. When 
compared to []l9| the only difference is in the mass term. 
Assuming m > in (|l|) leads to two types of single site 
interactions. The one with the solid bond is always posi- 
tive, and the one with the dotted bond is negative on filled 
even sites and empty odd sites. This extra negative sign 
must be included in Sign[n, b] along with sign £ that arises 
due to fermion hops, staggered fermion phase factors and 
anti-periodic spatial boundary conditions as discussed in 

Bonds connect lattice points into clusters. A flip of a 
cluster is defined as the change in the configuration of 
fermion occupation numbers at the sites belonging to the 
cluster, such that an occupied site is emptied and vice- 
versa. For a given model to be solvable using a meron 
cluster algorithm, the weights W[n, b] and signs Sign[n, 6] 
must satisfy three properties under cluster flips. 

1. The weight of the configuration W[n, b] must not 
change under the flip of any cluster. 

2. The change in Sign[n, b] due to a cluster flip must 
be independent of the state of other clusters. 



3. Starting from any configuration [n,b], it must be 
possible to reach a reference configuration [n re f,6] 
by flipping clusters such that Sign[n xe f,6] is 1. 

As has been discussed in [^5| , the change in the sign of the 
configuration due to a cluster flip depends on the topol- 
ogy of the cluster. If we refer to the clusters whose flip 
changes the sign of the configuration as merons, then it is 
easy to classify a configuration [n, b] based on the number 
of meron clusters it contains. Using the above three prop- 
erties it is then easy to check that the partition function 
gets contributions only from the zero meron sector, i.e., 

Z =^2Sign[n,b]W[n,b]. (7) 

[n,5] 

where Sign[n, b] = Sn.o and N denotes the number of 
merons in the configuration. For the present model it 
is also easy to show that the chiral condensate is given by 

(H>) = y^Y. Size(C mcron ) 6 NA W[n, b], (8) 

[n,b] 

which means that to measure the condensate one is inter- 
ested only in the zero and one meron number sectors. 

A typical sweep in a meron cluster algorithm consists of 
two steps exactly like any other cluster algorithm. Start- 
ing from a configuration the bonds are updated first based 
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on the weights W[n, b}. After each local bond update the 
meron number of the configuration can potentially change. 
In order not to generate more than the necessary number 
of merons every bond update is followed by a Metropolis 
decision. While measuring the chiral condensate for exam- 
ple one rejects all bond updates that generate more than 
one meron. This requires one to reanalyze the topology 
of clusters that are connected to the local bonds being 
updated. A major improvement in the implementation 
allows this to be done in at most log(Size(C)) steps [ po| . 
Once all the bonds are updated it is possible to update 
the occupation numbers by flipping each cluster with a 
probability half. This algorithm produces the configura- 
tion [n, b] with weight (Sn. + dN,i)W[n, b}. Thus, the 
condensate can be calculated using 
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(9) 



V( <5jv,o ) 

where (...) refers to a simple average over the generated 
configurations. 

4. Numerical Results 

The critical behavior of the model was studied through 
the measurement of the chiral condensate at several val- 
ues of the temperature around T c using the algorithm 
discussed above. For each temperature simulations were 
performed on different spatial lattice sizes and at various 
masses. Each simulation produced 10 6 configurations, ex- 
cept for the largest lattices of sizes ranging from 32 3 up 
to 48 3 which contained only 10 5 configurations. All runs 
included at least ten thousand thermalization sweeps in 
addition to the above measurements. The autocorrelation 
times typically ranged from two to five sweeps. However, 
errors were evaluated from fluctuations in the averages of 
data over blocks of 1000 configurations each. For future 
reference we give the values of the chiral condensate on a 
32 3 lattice obtained at m = O.OOf and various tempera- 
tures in the table below. 



T 




T 


W> 


1.0000 
1.0152 
1.0309 


0.270(20) 
0.228(16) 
0.205(12) 


1.0471 
1.0638 
1.0870 


0.1478(61) 
0.0566(17) 
0.0175(04) 



In order to confirm the spontaneous breaking of the Z2 
chiral symmetry the condensate needs to be evaluated in 
the infinite volume limit followed by the chiral limit. This 
can be done precisely by a finite size scaling analysis. In 
the broken phase the theory undergoes a first order phase 
transition as a function of the mass at m = where the 
condensate exhibits a jump. In a large but finite volume 
this discontinuity is smoothencd out to an analytic curve 
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Figure 3. Chiral condensate as a function of the mass at 
three different lattice sizes at T — 1.0471. The data from 
the simulations is plotted along with the function ( ^fl| ) for 
each volume. The same values of So and Xo are used for 
all curves shown. 



whose functional form is given by 
(vjjip) = E tanh(ml/S /T) + xo m , 
which is valid when 
m < S /xo • 



(10) 



(11) 



By fitting the available data at each simulation tempera- 
ture to the formula (10) it is possible to extract So, the 
desired limiting value of the condensate. The minimum 
volume necessary for the formula to work can be deter- 
mined by systematically removing the smallest volume 
data from the fit, as required to obtain a \ 2 /DOF of 
about one. Finally, it is important to check if So and xo 
obtained through the fit and the masses used are consis- 
tent with the condition (pi]). 

The above finite size scaling analysis works exception- 
ally well, as can be seen from figure which shows a plot 
of the condensate as a function of the mass at a fixed 
T = 1.0471 for three different lattice sizes. The solid 
lines represent the function ([!(]) for a fixed value of So 
and xo obtained from a single fit to all of the shown data 
and more. The x 2 /DOF for the fit is around 0.8. Since 
T = 1.0471 is very close to the critical temperature it was 
necessary to go to spatial volumes as large as 48 3 in order 
to determine So with reasonable precision. 

Interestingly, the above fitting procedure failed to yield 
an acceptable chi-squared for data above a certain tem- 
perature. This was essentially due to the fact that it was 
impossible to find a small enough mass to ensure that the 
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Figure 4. Chiral condensate versus temperature. The 
open squares are results for L = 32 and m — 0.001. The 
filled diamonds are the extrapolation to the infinite vol- 
ume limit and the chiral limit. The solid line is a plot of 
v4(1.0518 - T) 314 obtained from a fit of the extrapolated 
points. 
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Figure 5. Log-log plot of L d ~ Vm (i/jip) versus x = L Vm m 
which describes the scaling function given in (|l2|). The 
value of y rn used is obtained from a fit to the small x data 
as explained in the text. The solid lines are fits to (|l|) for 
small x and ( p^ ) for large x data. 



condition (|llj) was satisfied. Removing the larger mass 
data resulted in a smaller So which in turn lowered the 
bound in (|l"l|). This observation is consistent with a van- 
ishing condensate in the chiral limit for higher tempera- 
tures. 

The non-zero values of Eo extracted from the fits at 
different temperatures are shown in figure ^. The solid line 
represents a fit to the expected form A(T C — T) 13 using the 
data for temperatures between 1.0050 and 1.0471. The 
X 2 /DOF for the fit was 0.4. The critical exponent (3 was 
found to be 0.314(7) which is consistent with the measured 
value of 0.32648(18) in the Ising model §§] within two 
sigma. The fit also determines the critical temperature 
very accurately to be T c — 1.0518(3) which agrees with 
the previously obtained result |L9| within errors. The data 
for L = 32 and m — 0.001 is shown for comparison and to 
demonstrate the necessity of using the finite size scaling 
formula ([!(]) to extract the infinite volume chiral limit. 
Furthermore, the errors on Eo are typically much smaller 
than those for any single data point due to the constraints 
that arise from fitting over a wide range of masses and 
volumes. 

At T c the finite size scaling behavior of the condensate 
as a function of the mass is very different from (10) and 
can be shown to be of the form 



using renormalization group arguments that are valid close 
to a second order phase transition. For any Z2 phase 
transition the function f(x) is universal and depends only 
on the lattice geometry, boundary conditions etc., and not 
on the type of lattice, irrelevant operators and such. The 
exponent y m is hence universal. The behavior of f{x) is 
known in the two limits: 



/(*) 
fix) 



fo x , 

foo X~& , 



(13) 
(14) 



where 1/5 = (d - y m )/y m - 

Based on the estimate of T c obtained above, additional 
simulations at T = 1.0518 were performed to confirm this 
critical behavior as a function of the mass. Using the data 
with L > 8 and L Vm m < 0.1, the chiral condensate is fit 
to the form 



(15) 



based on the expectation from (|13|). Using this value of 
y m , we can calculate 



6 = 4.87(10) 



(16) 



(12) 



which compares favorably with the known value for the 
Ising model of 4.7893(22)_ |§. Using the fitted value of 
y m , the quantity L d ~ Vm {tpip) is plotted in figure || for all 
available data as a function of L Vm m. Clearly the data 
collapses onto a single curve within errors, as expected 
from ([12]) . To confirm this picture the data in figure || 



G 



with L > 12 and 10 < L Vm m < 90 is fit to the form (|14|). 
This yields the value 5 — 4.89(19) which agrees with the 
earlier determination and further demonstrates that the 
critical behavior in the present model is consistent with 
the 3-d Ising universality class. 

5. Conclusions 

The meron cluster algorithm has allowed a direct sim- 
ulation of a strongly interacting fermionic theory in the 
vicinity of a Z2 chiral phase transition on lattices with 
spatial volumes of up to 48 3 using common workstation 
computers. The scaling behavior near the critical point, 
as a function of both the temperature and mass, was de- 
termined within errors of about 2% and the results are 
consistent with a second order phase transition. Allowing 
for 1-2 sigma deviations the universality class of the transi- 
tion matches with that of the 3-d Ising model. This result 
strongly supports the scenario where a fermionic theory in 
d + 1 dimensions undergoes a dimensional reduction to be 
described by a d dimensional bosonic theory in the critical 
region. Comparing similar results with the hybrid Monte 
Carlo algorithm, the meron algorithm provides great im- 
provement over such methods. 
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